%------------------------------------------------------------------------%
%              Comunicaciones digitales y analógicas                     %
%                      Trabajo Práctico N 1                              %
%                                                                        %
%                       Eric Lionel Koplin                               %
%------------------------------------------------------------------------%
close all;clear all;clc;
%% Probabilidad de error de una M-PAM
%---------------------------------------------------------------Parámetros
% Parametros de ajuste de la simulacion
%
M = 4; % cantidad de simbolos
k = log2(M); % bits por simbolo
SNRbdB = 2:14; % SNR por bit en dB
nsimbs = 2e4 ; % cantidad de simbolos que para cada valor de SNRbdB en el Montercarlo

%----------------------------------- Conversiones e inicializaciones varias
%
SNRbveces = 10.^(SNRbdB/10); % SNR por bit en veces
SNRdB = SNRbdB + 10*log10(k);% SNR en dB
SNR = 10.^(SNRdB/10);% SNR en veces

Eavg = 1 ;
Pes = zeros(size(SNRbdB));% probabilidad de error por simbolo simulada
Peb = zeros(size(SNRbdB));% probabilidad de error de bit simulada

%--------------------------------------------------------P(error) teorica :
%
% Probabilidad de error de bit teorica de una 4−PAM
% Uso un puntero a una funcion cuyo parametro es SNR (en veces)
Pet = @(SNR) 2*(1-1/M) * qfunc(sqrt(3*SNR/(M^2-1))) ;

%-----------------------------------------------------Simulacion Montecarlo
%--------------------PAM
% modulador
modPAM = comm.PAMModulator(M,'SymbolMapping','Gray',...
                          'NormalizationMethod', 'Average Power',...
                          'AveragePower',Eavg);
%demodulador
demodPAM = comm.PAMDemodulator(M,'SymbolMapping','Gray',...
                              'NormalizationMethod' , 'Average Power' , ...
                              'AveragePower',Eavg);
%--------------------PSK
%Recordar la energia está normalizada, Es=Eavg=1
% modulador
modPSK = comm.PSKModulator(M,'SymbolMapping','Gray');% Eavg=1
%demodulador
demodPSK = comm.PSKDemodulator(M,'SymbolMapping','Gray');
%--------------------QAM
% modulador
modQAM = comm.RectangularQAMModulator(M,'SymbolMapping','Gray',...
                          'NormalizationMethod', 'Average Power',...
                          'AveragePower',Eavg);
%demodulador
demodQAM = comm.RectangularQAMDemodulator(M,'SymbolMapping','Gray',...
                              'NormalizationMethod' , 'Average Power' , ...
                              'AveragePower',Eavg); 
%--------------------FSK

%Recordar la energia está normalizada, Es=Eavg=1
% modulador
modFSK = comm.FSKModulator('ModulationOrder',M,'SymbolMapping','Gray','SamplesPerSymbol',1);% Eavg=1
%demodulador
demodFSK = comm.FSKDemodulator('ModulationOrder',M,'SymbolMapping','Gray','SamplesPerSymbol',1);                          

%--------------------AWGN Channel
hAWGN = comm.AWGNChannel('NoiseMethod','Signal to noise ratio (Eb/No)',...
    'EbNo',SNRbdB(1),'SignalPower',Eavg,'BitsPerSymbol',k);
%hAWGN = comm.AWGNChannel('NoiseMethod', ...
%                         'Signal to noise ratio (SNR)','SNR',0);
%--------------------Create an error rate calculator
hError = comm.ErrorRate;
%% recordar hacer sigma/2 en el teorico tambien
simbtx = randi(M,nsimbs*length(SNRbdB),1)-1;% simbolos a transmitir (en 1:M);
simbmod = step(modFSK,simbtx); % simbolos modulados 
%simbmod_ruido = step(hAWGN,reshape(simbmod,nsimbs,length(SNRbdB)));%en cada columba una canal distinto
sigma = sqrt(Eavg./SNR); % desvio de ruido
ruido = complex(randn(nsimbs,1),randn(nsimbs,1))*sigma;
simbmod_ruido = simbmod + reshape(ruido,nsimbs*length(SNRbdB),1);

simbrx= step(demodFSK,reshape(simbmod_ruido,nsimbs*length(SNRbdB),1)) ; % simbolos recibidos
%errores = step(hError, simbtx, simbrx);

%cada columna es un canal con distinta snr y cada elemento es una señal

Pes = ones(1,nsimbs)*reshape(simbtx ~= simbrx,nsimbs,length(SNRbdB))/nsimbs ;

% Ahora analizo los errores de bit (el Gray lo hizo el modulador )
%
bitstx = de2bi(simbtx,k) ; % convierto los simbolos transmitidos a bits
bitsrx = de2bi(simbrx,k) ; % convierto los simbolos recibidos a bits
% comparo los bits y cuento cuantos difieren
Peb = ones(1,nsimbs)*reshape(xor(bitstx,bitsrx)*ones(k,1),nsimbs,length(SNRbdB)) /k/nsimbs;

% Graficos
%
% Probabidad de error de simbolo
figure(3)
semilogy(SNRbdB,Pet(SNR),'-dr','linewidth',1.5);hold on
semilogy(SNRbdB,Pes,'-o','linewidth',1.5);hold off
grid on
box on
xlabel('SNR por bit (dB)') ;
ylabel('Pe')
legend(' Curva teorica','Montecarlo')

% Probabilidad de error por bit
figure(4)
semilogy(SNRbdB,Pet(SNR)/k,'-dr');hold on
semilogy(SNRbdB,Peb,'-o','linewidth',1.5);hold off
grid on
box on
xlabel('SNR por bit (dB)') ;
ylabel('Pb')
legend('Curva teorica','Montecarlo')
%{
%--------------------Ciclo Pe dada SNR                         
for n0 = 1:length(SNRbdB)
% Primero trabajo con simbolos y despues analizo por bit
%
sigma = sqrt(Eavg/SNR(n0)); % desvio de ruido
simbtx = randi(M,nsimbs,1)-1;% simbolos a transmitir (en 1:M);

% la funcion ”step” llama al objeto modulador o demodulador
simbmod = step(modPAM,simbtx); % simbolos modulados para su transmision ”OJO , son COMPLEJOS”
% Sumo el ruido (OJO sumo el ruido como si la constelacion tuviera parte compleja )
simbmod = simbmod + sigma * (randn(nsimbs,1) + 1i*randn(nsimbs,1));
% Demodulo
simbrx= step(demodPAM,simbmod) ; % simbolos recibidos
% Probabilidad de error de simbolo simulada
Pes(n0) = nnz(simbtx ~= simbrx )/nsimbs ;

% Ahora analizo los errores de bit (el Gray lo hizo el modulador )
%
bitstx = de2bi(simbtx) ; % convierto los simbolos transmitidos a bits
bitsrx = de2bi(simbrx) ; % convierto los simbolos recibidos a bits
% comparo los bits y cuento cuantos difieren
Peb(n0) = nnz(xor(bitstx,bitsrx)) /k/nsimbs ;
end


% Graficos
%
% Probabidad de error de simbolo
figure(1)
semilogy(SNRbdB,Pet(SNR),'-dr','linewidth',1.5);
hold on
semilogy(SNRbdB,Pes,'-o','linewidth',1.5) ;
grid on
box on
xlabel('SNR por bit (dB)') ;
ylabel('Pe')
legend(' Curva teorica','Montecarlo')

% Probabilidad de error por bit
figure(2)
semilogy(SNRbdB,Pet(SNR)/k,'-dr');
hold on
semilogy(SNRbdB,Peb,'-o','linewidth',1.5) ;
grid on
box on
xlabel('SNR por bit (dB)') ;
ylabel('Pb')
legend('Curva teorica','Montecarlo')

%%
%--------------------Ciclo Pe dada SNR                         
for n0 = 1:length(SNRbdB)
% Primero trabajo con simbolos y despues analizo por bit
%
sigma = sqrt(Eavg/SNR(n0)); % desvio de ruido
simbtx = randi(M,nsimbs,1)-1;% simbolos a transmitir (en 1:M);

% la funcion ”step” llama al objeto modulador o demodulador
simbmod = step(modFSK,simbtx); % simbolos modulados para su transmision ”OJO , son COMPLEJOS”
% Sumo el ruido (OJO sumo el ruido como si la constelacion tuviera parte compleja )
simbmod = simbmod + sigma * (randn(nsimbs,1) + 1i*randn(nsimbs,1));
% Demodulo
simbrx= step(demodFSK,simbmod) ; % simbolos recibidos
% Probabilidad de error de simbolo simulada
Pes(n0) = nnz(simbtx ~= simbrx )/nsimbs ;

% Ahora analizo los errores de bit (el Gray lo hizo el modulador )
%
bitstx = de2bi(simbtx) ; % convierto los simbolos transmitidos a bits
bitsrx = de2bi(simbrx) ; % convierto los simbolos recibidos a bits
% comparo los bits y cuento cuantos difieren
Peb(n0) = nnz(xor(bitstx,bitsrx)) /k/nsimbs ;
end



% Graficos
%
% Probabidad de error de simbolo
figure(3)
semilogy(SNRbdB,Pet(SNR),'-dr','linewidth',1.5);
hold on
semilogy(SNRbdB,Pes,'-o','linewidth',1.5) ;
grid on
box on
xlabel('SNR por bit (dB)') ;
ylabel('Pe')
legend(' Curva teorica','Montecarlo')

% Probabilidad de error por bit
figure(4)
semilogy(SNRbdB,Pet(SNR)/k,'-dr');
hold on
semilogy(SNRbdB,Peb,'-o','linewidth',1.5) ;
grid on
box on
xlabel('SNR por bit (dB)') ;
ylabel('Pb')
legend('Curva teorica','Montecarlo')

%%

% Example: 
    %   Modulate and demodulate a signal using 8-FSK modulation with a
    %   frequency separation of 100 Hz
 
     hMod = comm.FSKModulator(8, 100,'BitInput',false);
     hAWGN = comm.AWGNChannel('NoiseMethod', ...
                         'Signal to noise ratio (SNR)','SNR',0);
     hDemod = comm.FSKDemodulator(8, 100);
     %Create an error rate calculator
     hError = comm.ErrorRate;
     for counter = 1:100
         % Transmit a 50-symbol frame
         data = randi([0 hMod.ModulationOrder-1],50,1);
         modSignal = step(hMod, data);
         noisySignal = step(hAWGN, modSignal);
         receivedData = step(hDemod, noisySignal);
         errorStats = step(hError, data, receivedData);
     end
    fprintf('Error rate = %f\nNumber of errors = %d\n', ...
      errorStats(1), errorStats(2))
%}

